home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATH1 / GAUSID.PAS < prev    next >
Pascal/Delphi Source File  |  1985-04-03  |  3KB  |  179 lines
  1. program gausid;        { -> 129 }
  2. { pascal program to perform simultaneous solution }
  3. { by Gauss-Seidel }
  4. { procedure SEID is included }
  5.  
  6. const    maxr    = 8;
  7.     maxc    = 8;
  8.  
  9. type    ary    = array[1..maxr] of real;
  10.     arys    = array[1..maxc] of real;
  11.     ary2s    = array[1..maxr,1..maxc] of real;
  12.  
  13. var    y    : ary;
  14.     coef    : arys;
  15.     a    : ary2s;
  16.     n,m    : integer;
  17.     first,
  18.     error    : boolean;
  19.  
  20. external procedure cls;
  21.  
  22. procedure get_data
  23.     (var a    : ary2s;
  24.      var y    : ary;
  25.      var n,m: integer);
  26. { get values for n and arrays a,y }
  27.  
  28. var    i,j    : integer;
  29.  
  30. begin
  31.   writeln;
  32.   repeat
  33.     write('How many equations? ');
  34.     readln(n);
  35.     if first then first:=false else cls
  36.   until n<maxr;
  37.   m:=n;
  38.   if n>1 then
  39.     begin
  40.       for i:=1 to n do
  41.     begin
  42.       writeln('Equation',i:3);
  43.       for j:=1 to n do
  44.         begin
  45.           write(j:3,':');
  46.           read(a[i,j])
  47.         end;
  48.     write(' C:');
  49.     read(y[i]);
  50.     readln        { clear the line }
  51.       end;
  52.     writeln;
  53.     for i:=1 to n do
  54.       begin
  55.     for j:=1 to m do
  56.       write(a[i,j]:7:4,' ');
  57.     writeln(':',y[i]:7:4)
  58.       end;
  59.       writeln
  60.     end        { if n>1 }
  61.   else if n<0 then n:=-n;
  62.   m:=n
  63. end;        { procedure get_data }
  64.  
  65. procedure write_data;
  66. { print out the answers }
  67.  
  68. var    i : integer;
  69.  
  70. begin
  71.   for i:=1 to m do
  72.     write(coef[i]:9:5);
  73.   writeln
  74. end;        { write_data }
  75.  
  76. procedure seid
  77.     (a    : ary2s;
  78.     y    : ary;
  79.     var coef: arys;
  80.     ncol    : integer;
  81.        var error: boolean);
  82. { matrix solution by Gauss Seidel }
  83.  
  84. const    tol    = 1.0E-4;
  85.     max    = 100;
  86.  
  87. var    done    : boolean;
  88.        i,j,k,l,n: integer;
  89.  
  90.     nextc,hold,
  91.     sum,lambda,
  92.     ab,big    : real;
  93.  
  94. begin
  95.   repeat
  96.     write('Relaxation factor? ');
  97.     readln(lambda)
  98.   until (lambda<2) and (lambda>0.0);
  99.   error:=false;
  100.   n:=ncol;
  101.   for i:=1 to n-1 do
  102.     begin
  103.       big:=abs(a[i,i]);
  104.       l:=i;
  105.       for j:=i+1 to n do
  106.     begin
  107.     { search for largest element }
  108.       ab:=abs(a[j,i]);
  109.       if ab>big then
  110.         begin
  111.           big:=ab;
  112.           l:=j
  113.         end
  114.       end;        { j-loop }
  115.     if big=0.0 then error:=true
  116.     else
  117.       begin
  118.     if l<>i then
  119.       begin
  120.       { interchange rows to put }
  121.       { largest element on diagonal }
  122.       for j:=1 to n do
  123.         begin
  124.           hold:=a[l,j];
  125.           a[l,j]:=a[i,j];
  126.           a[i,j]:=hold
  127.         end;
  128.       hold:=y[l];
  129.       y[l]:=y[i];
  130.       y[i]:=hold
  131.     end    { if l<>i }
  132.       end    { if big }
  133.     end;    { i-loop }
  134.   if a[n,n]=0.0 then error:=true
  135.   else
  136.     begin
  137.       for i:=1 to n do
  138.     coef[i]:=0.0;        { initial guess }
  139.       i:=0;
  140.       repeat
  141.     i:=i+1;
  142.     done:=true;
  143.     for j:=1 to n do
  144.       begin
  145.         sum:=y[j];
  146.         for k:=1 to n do
  147.           if j<>k then
  148.         sum:=sum-a[j,k]*coef[k];
  149.         nextc:=sum/a[j,j];
  150.         if abs(nextc-coef[j])>tol then
  151.           begin
  152.         done:=false;
  153.         if nextc*coef[j]<0.0 then
  154.           nextc:=(coef[j]+nextc)*0.5
  155.           end;
  156.         coef[j]:=lambda*nextc+(1.0-lambda)*coef[j];
  157.         writeln(i:4,',coef(',j,')=',coef[j])
  158.     end        { j-loop }
  159.       until done or (i>max)
  160.   end;    { if a[n,n]=0 }
  161.   if i>max then error:=true;
  162.   if error then writeln('ERROR: Matrix is singular')
  163. end;        { SEID }
  164.  
  165. begin        { MAIN program }
  166.   first:=true;
  167.   cls;
  168.   writeln;
  169.   writeln('Simultaneous solution by Gauss-Seidel');
  170.   repeat
  171.     get_data(a,y,n,m);
  172.     if n>1 then
  173.       begin
  174.     seid(a,y,coef,n,error);
  175.     if not error then write_data
  176.       end
  177.   until n<2
  178. end.
  179.